#!/usr/bin/perl

use strict;
use PGPLOT;
use POSIX;

# This script makes a plot to enable a comparison between effectiveness of
# on-the-fly (OTF) and regular mosaicing strategies.

# Jamie Stevens, 2010-12-02

# go through the arguments
my %arguments;
$arguments{"plot_device"}="/xs";
$arguments{"frequency"}=1400.0;
$arguments{"min_cycles"}=2;
for (my $i=0;$i<=$#ARGV;$i++){
    if ($ARGV[$i] eq "--plot"){
	$i++;
	$arguments{"plot_device"}=$ARGV[$i];
    } elsif ($ARGV[$i] eq "--freq"){
	$i++;
	$arguments{"frequency"}=$ARGV[$i];
    } elsif ($ARGV[$i] eq "--spec-size"){
	$i++;
	$arguments{"specific_size"}=$ARGV[$i];
    } elsif ($ARGV[$i] eq "--min-cycles"){
	$i++;
	$arguments{"min_cycles"}=$ARGV[$i];
    }
}

# get the viewport size
#pgopen("mosaicing_comparison.ps/cps");
pgopen($arguments{"plot_device"});
my ($vp_x1,$vp_x2,$vp_y1,$vp_y2);
pgqvp(0,$vp_x1,$vp_x2,$vp_y1,$vp_y2);

# the range for numbers of individual pointings
my ($min_n_pointings,$max_n_pointings,$incr_n_pointings)=(100,400,10);
my @n_pointings;
for (my $i=$min_n_pointings;$i<=$max_n_pointings;$i+=$incr_n_pointings){
    push @n_pointings,$i;
}

# the range for correlator cycle times (seconds)
#my ($min_cycle,$max_cycle,$incr_cycle)=(2,10,4);
my @cycles=(2,10);

# the range of number of cycles spent at each pointing in normal mosaicing
my ($min_n_cycles,$max_n_cycles,$incr_n_cycles)=($arguments{"min_cycles"},4,2);

# the number of hours in total
my $total_hours=12.0;

# the frequency to compute for (in MHz)
my $frequency=$arguments{"frequency"};

# the speed of light (m/s)
my $c=299792458.0;

# the diameter of the dish (m)
my $dish_diameter=22.0;

# the radian to degree conversion
my $rad_to_deg=180.0/3.141592654;

# the acceleration rate of the dish (deg/min/min)
my $accel_rate=800.0;

# start computations
# how far apart are the beam centres?
my $pointing_separation=(0.5*($c/($frequency*1e6))/$dish_diameter)*
    $rad_to_deg; # in degrees

# how long does it take to slew between pointing centres (normal
# mosaicing)?
my $nm_slew_time=2.0*sqrt($pointing_separation/$accel_rate); # in minutes
print "slewing time between pointing centres: ".$nm_slew_time." (minutes)\n";

# start the plot
my ($min_time,$max_time)=(0,log(10)/log(10));
# x: pointings, y: time (min)
&panel(0);
pglab("","","Mosaicing at ".$frequency." MHz");
&panel(1);
#pgswin($min_n_pointings,$max_n_pointings,$min_time,$max_time);
pgbox("BCNTS",0.0,0,"BCLNTS",0.0,0);
pglab("Number of pointings","Time per pointing after ".$total_hours.
      " hours (minutes)","");
my ($min_fieldtime,$max_fieldtime)=(log(2)/log(10),log(400)/log(10));
&panel(2);
#pgswin($min_n_pointings,$max_n_pointings,$min_visits,$max_visits);
pgbox("BCTS",0.0,0,"BCLMTS",0.0,0);
pgmtxt("R",2.5,0.5,0.5,"Field survey time (minutes)");
#pglab("","Number of visits per pointing","");
#pgswin($min_n_pointings,$max_n_pointings,$min_time,$max_time);
&panel(1);

my $x_colour_label_left=$min_n_pointings+($max_n_pointings-$min_n_pointings)/10.0;
my $y_colour_label_top=$max_time-($max_time-$min_time)/10.0;
my $y_colour_label_incr=-1*($max_time-$min_time)/10.0;

my $x_line_label_left=$max_n_pointings-($max_n_pointings-$min_n_pointings)/6.0;
my $y_line_label_top=$y_colour_label_top;
my $y_line_label_incr=$y_colour_label_incr;
my $x_line_right_coord=$x_line_label_left-
    ($max_n_pointings-$min_n_pointings)/20.0;
my $x_line_left_coord=$x_line_right_coord-
    ($max_n_pointings-$min_n_pointings)/10.0;

my $curr_colour;
my $curr_linestyle=1;
foreach my $j (@cycles){
    # go through each cycle time
    pgsci(1);
    pgsls($curr_linestyle);
    # make the linestyle label
    my @x_line_locs=($x_line_left_coord,$x_line_right_coord);
    my @y_line_locs=($y_line_label_top+($curr_linestyle-1)*$y_line_label_incr,
		     $y_line_label_top+($curr_linestyle-1)*$y_line_label_incr);
    pgline(2,\@x_line_locs,\@y_line_locs);
    pgtext($x_line_label_left,$y_line_locs[0],$j." seconds");
    $curr_colour=1;
    for (my $k=$min_n_cycles;$k<=$max_n_cycles;$k+=$incr_n_cycles){
	my $spec_done=0;
	if (!$arguments{"specific_size"}){
	    $spec_done=-1;
	}
	pgsci($curr_colour);
	# make the colour label
	pgtext($x_colour_label_left,$y_colour_label_top+
	       ($curr_colour-1)*$y_colour_label_incr,$k." cycles");
	# go through the number of cycles on each pointing
	# (normal mosaicing)
	my @this_n_pointings;
	my @this_pointing_int_time;
	my @this_fieldtime;
	for (my $i=0;$i<=$#n_pointings;$i++){
	    # cycle through the number of pointings
	    my $spec_region_side_length;
	    my $spec_n_side_pointings;
	    my $spec_n_pointings;
	    if ($spec_done==0){
		$spec_region_side_length=sqrt($arguments{"specific_size"});	
		$spec_n_side_pointings=$spec_region_side_length/$pointing_separation;
		$spec_n_pointings=$spec_n_side_pointings*$spec_n_side_pointings;
	    }
	    # how long are we on each pointing
	    my $each_pointing_time=($j*$k)/60.0; # in minutes
	    # how much time for each total field?
	    my $total_field_time=$each_pointing_time*$n_pointings[$i];
	    my $spec_total_field_time;
	    if ($spec_done==0){
		$spec_total_field_time=$each_pointing_time*$spec_n_pointings;
	    }
	    # how much total time slewing?
	    my $total_slewing_time=($n_pointings[$i]-1)*$nm_slew_time;
	    my $spec_total_slewing_time;
	    if ($spec_done==0){
		$spec_total_slewing_time=($spec_n_pointings-1)*$nm_slew_time;
	    }
	    # how much integration time per field then?
	    my $total_field_integration_time=$total_field_time-
		$total_slewing_time;
	    my $spec_total_field_integration_time;
	    if ($spec_done==0){
		$spec_total_field_integration_time=$spec_total_field_time-
		    $spec_total_slewing_time;
	    }
	    # how many fields can we do in the total time?
	    my $total_n_fields=($total_hours*60.0)/$total_field_time;
	    my $spec_total_n_fields;
	    if ($spec_done==0){
		$spec_total_n_fields=($total_hours*60.0)/$spec_total_field_time;
	    }
	    # what integration time does that correspond to?
	    my $total_integration_time=$total_field_integration_time*
		$total_n_fields;
	    my $spec_total_integration_time;
	    if ($spec_done==0){
		$spec_total_integration_time=$spec_total_field_integration_time*
		    $spec_total_n_fields;
	    }
	    # what integration time per pointing does that correspond to?
	    my $total_pointing_integration_time=$total_integration_time/
		$n_pointings[$i];
	    my $spec_total_pointing_integration_time;
	    if ($spec_done==0){
		$spec_total_pointing_integration_time=$spec_total_integration_time/
		    $spec_n_pointings;
		print "\n for ".$arguments{"specific_size"}." square degrees (".
		    $spec_n_pointings." pointings), ".$j."s cycle time, ".$k." cycles:\n";
		print "  pointing integration time: ".$spec_total_pointing_integration_time." min/12 hours \n";
		print "  field survey time: ".$spec_total_field_time." min \n";
		$spec_done=1;
	    }
	    if ($total_pointing_integration_time>0){
		push @this_pointing_int_time,
		log($total_pointing_integration_time)/log(10.0);
		push @this_n_pointings,$n_pointings[$i];
		push @this_fieldtime,log($total_field_time)/log(10);
	    }
	}
	if ($#this_n_pointings>=0){
	    &panel(1);
#	    pgswin($min_n_pointings,$max_n_pointings,$min_time,$max_time);
	    pgline($#this_n_pointings+1,\@this_n_pointings,
		   \@this_pointing_int_time);
#	    pgswin($min_n_pointings,$max_n_pointings,$min_visits,$max_visits);
	    &panel(2);
	    pgline($#this_n_pointings+1,\@this_n_pointings,
		   \@this_fieldtime);
#	    pgswin($min_n_pointings,$max_n_pointings,$min_time,$max_time);
	    &panel(1);
	}
	$curr_colour++;
    }
    $curr_linestyle++;
}

my $x_colour_label_left=$min_n_pointings+3*($max_n_pointings-$min_n_pointings)/10.0;

$curr_linestyle=1;
foreach my $j (@cycles){
    pgsci($curr_colour);
    pgsls($curr_linestyle);
    my @this_f_n_pointings;
    my @this_f_pointing_int_time;
    my @this_f_fieldtime;
    # what is the slew rate for OTF mosaicing?
    my $slew_rate=$pointing_separation/($j/60); # degrees/min
    # the slew rate label
    my $slew_rate_label=sprintf "%.2f deg/min",$slew_rate;
    pgtext($x_colour_label_left,$y_colour_label_top+
	   ($curr_colour-3)*$y_colour_label_incr,
	   $slew_rate_label);
    my $spec_done=0;
    if (!$arguments{"specific_size"}){
	$spec_done=-1;
    }
    for (my $i=0;$i<=$#n_pointings;$i++){
	# how many pointings per line (assuming square region)?
	my $pointings_per_line=sqrt($n_pointings[$i]);
	my $spec_line_distance;
	my $spec_pointings_per_line;
	my $spec_n_pointings;
	if ($spec_done==0){
	    $spec_line_distance=sqrt($arguments{"specific_size"});	
	    $spec_pointings_per_line=$spec_line_distance/$pointing_separation;
	    $spec_n_pointings=$spec_pointings_per_line*$spec_pointings_per_line;
	}
	# how far along each line?
	my $line_distance=$pointings_per_line*$pointing_separation;
	my $closeness=($line_distance-floor($line_distance))/$line_distance;
#	print $line_distance." - ".floor($line_distance)." (".
#	    $closeness.")\n";
	my $tolerance=0.02;
	if (($closeness>=(0-$tolerance))&&
	    ($closeness<=($tolerance))){
	    &panel(2);
	    my @x_line_locs=($n_pointings[$i],$n_pointings[$i]);
	    my @y_line_locs=($min_fieldtime,$max_fieldtime);
	    pgsls(4);
	    pgsci(1);
	    pgline(2,\@x_line_locs,\@y_line_locs);
	    my $y_labloc=$min_fieldtime+($max_fieldtime-$min_fieldtime)/10.0;
	    pgtext($n_pointings[$i]+2,$y_labloc,floor($line_distance).'\uo\d x '.
		   floor($line_distance).'\uo\d');
	}
	# how long to do each line?
	my $line_integration_time=$line_distance/$slew_rate;
	my $spec_line_integration_time;
	if ($spec_done==0){
	    $spec_line_integration_time=$spec_line_distance/$slew_rate;
	}
	# assume a 2 second ramp time
	my $line_total_time=$line_integration_time+(2/60);
	my $spec_line_total_time;
	if ($spec_done==0){
	    $spec_line_total_time=$spec_line_integration_time+(2/60);
	}
	# multiply by the number of lines to do
	my $total_line_field_time=$line_total_time*$pointings_per_line;
	my $total_line_integration_time=$line_integration_time*
	    $pointings_per_line;
	my ($spec_total_line_field_time,$spec_total_line_integration_time);
	if ($spec_done==0){
	    $spec_total_line_field_time=$spec_line_total_time*
		$spec_pointings_per_line;
	    $spec_total_line_integration_time=$spec_line_integration_time*
		$spec_pointings_per_line;
	}
	# how many fields can we observe per total time?
	my $total_n_line_fields=($total_hours*60.0)/$total_line_field_time;
	my $spec_total_n_line_fields;
	if ($spec_done==0){
	    $spec_total_n_line_fields=($total_hours*60.0)/$spec_total_line_field_time;
	}
	# how much integration time does this correspond to?
	my $total_integration_time=$total_line_integration_time*
	    $total_n_line_fields;
	my $spec_total_integration_time;
	if ($spec_done==0){
	    $spec_total_integration_time=$spec_total_line_integration_time*
		$spec_total_n_line_fields;
	}
	# what integration time per pointing does that correspond to?
	my $total_pointing_integration_time=$total_integration_time/
	    $n_pointings[$i];
	my $spec_total_pointing_integration_time;
	if ($spec_done==0){
	    $spec_total_pointing_integration_time=$spec_total_integration_time/
		$spec_n_pointings;
	    print "\n for ".$arguments{"specific_size"}." square degrees (".
		$spec_n_pointings." pointings), ".$slew_rate_label.":\n";
	    print "  pointing integration time: ".$spec_total_pointing_integration_time." min/12 hours\n";
	    print "  field survey time: ".$spec_total_line_field_time." min\n";
	    $spec_done=1;
	}
	if ($total_pointing_integration_time>0){
	    push @this_f_pointing_int_time,
	    log($total_pointing_integration_time)/log(10.0);
	    push @this_f_n_pointings,$n_pointings[$i];
	    push @this_f_fieldtime,log($total_line_field_time)/log(10);
	}
    }
    if ($#this_f_n_pointings>=0){
	pgsci($curr_colour);
	pgsls($curr_linestyle);
	&panel(1);
	pgline($#this_f_n_pointings+1,\@this_f_n_pointings,
	       \@this_f_pointing_int_time);
	&panel(2);
	pgline($#this_f_n_pointings+1,\@this_f_n_pointings,
	       \@this_f_fieldtime);
	&panel(1);
    }
    $curr_colour++;
    $curr_linestyle++;
}

sub panel {
    my ($panel_number)=@_;

    if ($panel_number==1){
	# the bottom panel
	my $p_vp_y1=$vp_y1;
	my $p_vp_y2=$vp_y1+($vp_y2-$vp_y1)/2.0;
	pgsvp($vp_x1,$vp_x2,$p_vp_y1,$p_vp_y2);
	pgswin($min_n_pointings,$max_n_pointings,$min_time,$max_time);
    } elsif ($panel_number==0){
	# the whole panel
	pgsvp($vp_x1,$vp_x2,$vp_y1,$vp_y2);
    } elsif ($panel_number==2){
	# the top panel
	my $p_vp_y2=$vp_y2;
	my $p_vp_y1=$vp_y1+($vp_y2-$vp_y1)/2.0;
	pgsvp($vp_x1,$vp_x2,$p_vp_y1,$p_vp_y2);
	pgswin($min_n_pointings,$max_n_pointings,$min_fieldtime,$max_fieldtime);
    }
}
